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Abstract 

In order to clarify the statistical features of complex networks, the 
spectral density of adjacency matrices has often been investigated. 
Adopting a static model introduced by Goh, Kahng and Kim, we 
analyse the spectral density of complex scale free networks. For that 
purpose, we utilize the replica method and effective medium approxi- 
mation (EM A) in statistical mechanics. As a result, we identify a new 
integral equation which determines the asymptotic spectral density of 
scale free networks with a finite mean degree p. In the limit p — > do, 
known asymptotic formulae are rederived. Moreover, the 1/p correc- 
tions to known results are analytically calculated by a perturbative 
method. 
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1 Introduction 



The theory of complex networks, which originates in Erdos and Renyi's study 
on random graphs [1] , has dramatically developed in the last decade. Complex 
networks in biological and social systems were extensively investigated and 
universal features were elucidated. Moreover, as Erdos and Renyi's random 
graphs could not reproduce some of those features, new theoretical models 
of complex networks were introduced [2j [3]. 

One of the newly discovered features of complex networks is the scale 
free property, which means that the degree (the number of the nodes directly 
connected to each node) distribution obeys a power law. The scale free prop- 
erty was originally explained by Barabasi and Albert in terms of a growing 
network model[3]. Later Goh, Kahng and Kim showed that a static network 
model with the scale free property could also be constructed [1]. 



The connection patterns of the nodes in a complex network are repre- 
sented by the adjacency matrix. The spectral (eigenvalue) density of adja- 
cency matrices is used to measure the statistical features of large complex 
networks. For example, as explained in the following, a power law behaviour 
can be observed. Let us denote by N the size of the adjacency matrix. In 
the limit iV — > oo, we obtain an asymptotic form of the spectral density. 
When the mean value p of the degree becomes large, Erdos and Renyi's ran- 
dom graphs tend to have the asymptotic spectral density with the semicircle 
shape, as expected from a classical result of random matrix theory[5]. On the 
other hand, the tail of the asymptotic spectral density for scale free networks 
has a power law behaviour. This power law behaviour was analytically con- 
firmed by Dorogovtsev et al. for tree-like scale free graphsJHl [7J. Rodgers et 
al. also found a similar power law for Goh, Kahng and Kim's static model[5]. 

Rodgers et al. analysed the spectral density of the static model in the 
limit p — > oo [SI [9] . In their analysis, the replica method in statistical physics 
was employed. A similar usage of the replica method can be traced back to 
Rodgers and Bray's study of sparse random matrices |10j. Recently Semerjian 
and Cugliandolo[TT] proposed a scheme called the effective medium approx- 
imation (EMA) and elegantly rederived the results in [TO]. This scheme was 
also applied to sparse sample covariance matrices and analogous results were 
obtained [T2]. 

In this paper, we extend the works in [8] and pj] by treating complex 
networks with a finite mean degree p. For that purpose, we make use of the 
replica method and EMA for the analysis of Goh, Kahng and Kim's static 
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model. Consequently, in the asymptotic limit iV — > oo, we obtain a new 
integral equation which determines the asymptotic spectral density with a 
finite p. We show that known results originally found in [8] are rederived in 
the limit p — > oo. Moreover, using the new integral equation, we evaluate 
the 1/p corrections to the known results. 

This paper is organised as follows. In §2, the replica method is intro- 
duced for Goh, Kahng and Kim's static model of scale free networks. In §3, 
using EMA, we derive a nonlinear integral equation which determines the 
asymptotic spectral density with a mean degree p. In the limit p — > oo, the 
results in [8] are rederived. In §4, a perturbative method is developed to cal- 
culate the 1/p corrections. In §5 and §6, perturbation terms are analytically 
evaluated around the band center and in the tail region. 

2 Replica Method 

Goh, Kahng and Kim's static model of complex networks is constructed 
by the following procedurejl]. Let us suppose that there are N nodes and 
consider the asymptotic limit N —>■ oo. Each node j is assigned a probability 
Pj, which is normalised as YljLiPj = 1- in order to realize a power- law 
degree distribution with an exponent A = (1 + a) /a, we assume that 

Pj= ^i - - (l-a^r, 0<a<l. (2.1) 

In one step of the procedure, nodes j and I are chosen with probabilities Pj 
and Pi, respectively. Then the nodes j and I are connected by an edge, if j 
does not coincide with / and the nodes are not already connected. Repeating 
such a step pN/2 times, one obtains a network with a mean degree p. The 
probability that a pair of nodes j and I are connected by an edge is 

f jt = 1 - (1 - 2P j P l ) pN / 2 ~ 1 - e -P NP i p '. (2.2) 

In this paper, we are interested in the spectral density of the adjacency 
matrix of the network. The adjacency matrix J is symmetric (Jji = Jij) and 
the elements are independently distributed. The p. d.f. (probability distribu- 
tion function) of each element is given by 

V jl {J jl ) = {l-f jl )6{J jl ) + f jl 6{J jl -l), j<l. (2.3) 
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Denoting the average over the above probability distribution by brackets (•), 
we define the spectral density of J as 

P(/^) = (^X>(/^)}> (2-4) 

where /ii, /i2, • ■ • , Hn are the eigenvalues of J and S(x) is Dirac's delta func- 
tion. 

In order to analyse the asymptotic behaviour of the spectral density in 
the limit N — > oo, the partition function 

r N ( i N i N \ 

Z = JT[ dfy exp Ux ]T % - - £ Jjrfjh (2.5) 

is useful. In terms of Z, the spectral density can be rewritten as 
fa) = ^ImlV^J-^ + ie)/}- 1 ; 

= iv^ Im |l (lnZ(/i + ^ )) - (2 ' 6) 

Here / is an identity matrix and e is a positive infinitesimal number. More- 
over, the relation 

lim ln(^) = 

n->0 n 

can be employed to find 

"<" ) = I i,n o7vL Im | ln<{Z( ' l)} " > ' (2 ' 8) 

Hence we need to evaluate the average of the n-th power of Z. In terms 
of the replica variables 

6- = (4MV-,$°) (2-9) 
and the corresponding measures 

(ty,- = d^d^f-.-d^, (2.10) 

it can be written as 

(Z n ) = J ftdfj exp ^£$) ( ex P ^E^'-^)> (2-11) 
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where the average is evaluated as 

( ex p (-^e Jufi ■ A } = n (/ d ^ Wite-"** 1 *] 

N 

= n{ i +^( e "^'^- 1 )}- ( 2 - 12 ) 

Then one obtains an asymptotic estimate in the limit N — > oo as 

ln ( ex p E ^ • &\ ) ~ p^E ( e ~^" fi - x ) • ( 2 - 13 ) 

It was shown in [13] that the remainder term of this asymptotic estimate 
was 0(1) for < a < 1/2, 0((lniV) 2 ) for a = 1/2 and 0(N 2 ~^ IniV) for 
1/2 < a < 1. Moreover we define 

Cj$) = 8($-$j), (2.14) 

so that 

^exp ^E J H$i ■ 

N 



eip yE Pi* / #/d^ ~cS)~ciW (e-^ - l) | . (2.15) 



ii 

Then, introducing auxiliary functions Cj (</>), we find 

JV ( a N 



(n ~ /ii<%«pUa*e/« 

3=1 I 3=1 

X ex P |^E P 3^ Jd$Jdj ~ Cj {$)cS) (e- 

V JV V 

n^y n^wnn*few-c 3 w)e si+S2 , (2.16) 



where 



3=1 3=1 3=1 



A? 



5i = ^E/# c i(^ (2-17) 
z 3=1 7 



and 



S 2 = ^E Prf / d0 j d^ Cj {$) *$) (e-*W _ i) . 
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(2.18) 



Here the functional integration is taken over Cj(<f) satisfying 



d(p Cj(4>) = 1. 



(2.19) 



Let us note 



N N 

3=1 3=1 <a 



N r N 

3=1 3=1 



JV 



Y\_ T>aj{4>) exp 

3=1 



2 7 r^/d0a i (0){c i (0)-c i (0)} 

3=1 

3=1 L J 



2.20) 



where 



— In y d0j exp j— 27U y d</> aj((f))cj((f)) 

— In J d(f)j exp |— 27riaj(0j)| . 



(2.21) 



In the limit iV — ► oo, the functional integral ( 12.201) is dominated by the 
contribution from the neighbourhood of the extremum satisfying 



5 



Saj((j)) 
so that 



We thus obtain 

N N 



2iri [ d$ aj ($)cj($) - FA = 2m Cj ($) - 2me- 2ma ^ )+F] = 0, 



d(f) Cj(4>) h\Cj{4>) = 2ni / d<fi aj(<j>)cj(<j>) — Fj. 



[2.22) 
(2.23) 



N N j N I 

/ n <% nn%(^ - ~ e *p -z / # , (2.24) 

3=1 3=1 (A I 3=1 J 
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from which it follows that 

N 



Z n ) ~ J f[ Vc0) e So+Sl+S2 (2.25) 



with 

^o = -E /d0 9(0)1119(0). (2.26) 

3 Integral Equation 

In the limit iV — > oo, the dominant contribution to the functional integral 
(I2.25P again comes from the neighbourhood of the extremum. In order to 
work out the extremum, we replace Cj(0) with a Gaussian ansatz 

Then the extremum condition can be written as 

-£-(S + S 1 + S 2 )=0. (3.2) 
da k 

This variational scheme based on the Gaussian ansatz is called the effective 
medium approximation (EMA)[llJ. 
Using the identity 

d0 2 Cj(0) = nitjj, (3.3) 



we can readily find 



„ n . , n . , Nn n fin 

Z 3=1 Z Z j=l 



so that 

dS n dS\ fin 
da k 2ak da k 2 
It is straightforward to obtain 



(3.5) 



Np 
Np 



N 



1 



. x .. {2mai) 

N 



n/2 



cxp 



1 " 



2%ai 



PjPi 



(3.6) 



Therefore 



dS 2 Npn p * 



do, 



Then the extremum condition (13.21) yields 



JV 



„• 1 1 — cr,-crfc 



in the limit n — > 0. Putting (12. ip into (13.81) . we find 



(3.7) 



(3.8) 



= p p(l — a) 



kY a 1 " 



N N 



i=i 



(3.9) 



In the limit iV — > oo, it follows that 



cr x x 



px a -p(l-a) 2 [\ 
Jo 1 



a(x)a(y) 



dy 



(3.10) 



with x = k/N. The function cr(x) is determined by this integral equation. 
The spectral density can be evaluated as 

M = &i^ Im | ln<{Z( " )r> 



2 d 
lim — — Im— (S + S 1 + S 2 ) 

n^o JSn-n op 



A' 



Nix 



(3.11) 



i=i 



so that 



1 f 1 
pip) = Im / a(x)dx 

71 JO 



(3.12) 
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in the limit iV — > oo. The asymptotic spectral density is thus written in 
terms of the solution of the integral equation (I3.10p . 

Using (12.11) and (|2.2p . we generated adjacency matrices J with N = 1000 
and evaluated the spectral density (NUMERICAL) by numerical diagonal- 
ization. On the other hand, the EMA spectral density (EMA) was calculated 
by solving (13.91) with N = 1000 by numerical iteration and by putting the 
solution into (13. lip . We show the both results for p = 10 and a = 0.5 in 
Figure 1. The agreement seems good enough to demonstrate the validity of 
EMA. 
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Figure 1: Comparison of a numerically generated spectral density (average 
over 100 samples, = 1000) and an EMA solution (p = 10, a = 0.5). 
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4 Perturbation Theory 

In order to analytically evaluate the EMA spectral density, we start from 
the integral equation (13.101) and develop a perturbation theory for large p. 
Introducing scalings 

a(x) -> —s(x), fi -> y/pm, (4.1) 
VP 



we obtain 



s{x)x~ a = - ^— . (4.2) 

1 s(x)s(y) 

V 



Let us expand s{x) in terms of 1/p as 



s(x) = s (x) + -s 1 (x) + O[\), (4.3) 
V \P 



so that the unperturbed term Sq(x) satisfies 



s (x)x- a = — — . (4.4) 

mx a - (1 - af / y~ a s (y)dy 



This integral equation for the unperturbed term was originally found in [8] . 
Then, expanding the RHS of (14.21) in terms of 1/p, we find the relation 

^s (x) + ^ s i( x )\ x ~ a 

(1 - a) 2 r 1 r i 

~ s (x)x- a + [ - L [sQ{x) y x -^ / \ Sl (y) + s (x){s (y)} 2 }y~ a dy 

p Jo 1 J 

(4.5) 

in the limit p — > oo. Therefore, the perturbation term si(x) can be evaluated 

as 

Sl (x)x- a = (l-a) 2 {s (x)} 2 x- 2a f 1 Sl (y)y- a dy 

Jo 

+ (l-a) 2 {s (x)} 3 x- 2 « A So (y)}V a ck/, (4.6) 

Jo 
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so that 



Here 



and 



[ s 1 (x)dx= {l-afh{D + Ji). (4.7) 

JO 

J x = / 1 {s (a;)} 2 x- a dx, J x = / 1 {s (x)} 3 a;- Q dx (4.8) 
Jo Jo 



D = f 1 Sl (x)x~ a dx. (4.9) 
Jo 

Integrating the both sides of (14.61) over x from to 1, we find 

D = (l-a) 2 D ^{s^Yx'^dx + il-afh [' {s (x)} 3 x~ 2a dx 
Jo Jo 

= (l-a) 2 DI 2 + (l-a) 2 I 1 J 2 , (4.10) 

where 

h = f 1 {s {x)} 2 x~ 2a dx, J 2 = f '{s {x)} 3 x- 2a dx. (4.11) 
Jo Jo 

Therefore D can be evaluated as 

D= ilz^Ml (412) 

so that 

/; SlM dx = (l- a fA{ ^-^ 2 +Jl }. (4.13) 

Let us summarise the scheme of the perturbation. First we evaluate Sq(x) 
by solving the integral equation (14.41) . Next we calculate 

Im / s (x)dx (4.14) 

JO 

and the other integrals Ii, Ji, I 2 and J 2 . Then, using the relation (14.131) . we 
are able to compute 

Im / si{x)dx. (4-15) 
Jo 

In terms of the integrals (14.141) and (I4.15p . the asymptotic spectral density 
is expressed as 

p(fj) ~ — I Im / s (x)dx + -lmf si(x)6x + O ( 4r ) \ (4.16) 

71 yP [Jo p Jo \P J J 
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in the limit N — > oo. 

In the following sections, we evaluate the small and large m expansions of 
the integrals (14.141) and (I4.15p . so that (I4.16P gives the asymptotic behaviour 
of the EMA spectral density around the band center and in the tail region. 

5 Band Center 

Let us evaluate the EMA spectral density around the band center, by calcu- 
lating the small m expansions of the integrals (I4.14p and (I4.15p . Noting the 
dependence of s (x) on m, we expand it as 

s (x) = g (x) + mgi(x) + m 2 g 2 (x) + O (m 3 ) . (5.1) 

In order to compute the central value go(x), we put m = in (14. 4p and 
find 

g (x)x~ a = -I . (5.2) 

(1 - a) 2 / y- a g (y)dy 



o 

Integrating the both sides in x from to 1, we obtain 

i _ 1 

y a go(y)dy = -h , (5.3) 

o I — a 

so that 

go(x) = -t-^—x a . (5.4) 
1 — a 

Let us next evaluate gi(x) and g2{x). We expand (14.41) in terms of m and 
find 



[g (x) + mg^x) + m 2 g 2 (x)}x' 



a 



g {x)x 



-a 



1 - mg (x)x~ a |x Q - (1 - a) 2 £ gi(y)y- a dy} 



+ m 2 g (x)x a (l-a) 2 f g 2 {y)y a dy 

Jo 

+ m 2 {g (x)} 2 x- 2a [x a - (1 - a) 2 J* 9l (y)y- a dy 



(5.5) 
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in the limit m — > 0. The linear term in m reads 

^(x)^ = -{g ( x )} 2 x- 2a ^x a - (1 - a) 2 j\ 1 (y)y- a dy^ 

= tt-^ ~ f 9x{y)y- a &y- (5.6) 

(1 — a) z Jo 

We again integrate the both sides in x from to 1 and find 

/o lgl(y)y " ady= 2(l-a)'(l + a ) ' (5 - 7) 

so that 

9l{x)X ~ a = JTTaf ~ 2(1 -<*)*(! + a)' ^ 
The quadratic term in m is similarly extracted from (15.51) as 

g 2 (x)x- a = {g (x)} 2 x' 2a (l-a) 2 f 1 g 2 (y)y- a dy 



o 



{g (x)} 3 x- 3a {x a -(l-a) 2 f 1 gi (y)y- a dy }2 



o 



lo * y * Ja * (1-a) 3 [ 2(1 + a) 
Integrating the both sides in x from to 1 yields 

3 1 



2a 4(1 + a) 



(5.10) 



from which it follows that 



Thus we obtain the small m expansion of Sq(x): 

3 TYl f 1 I 

s (x)x~ a = —i- V — — < x° 



I- a (1 -a) 2 { 2(1 + a) J 

m 2 { 2n 2x a 5 1 1 ] / „x 
2x 2a + -- - } + O (m 3 ) , 



2(1 - a) 3 [ 1 + a 4(1 + a) 2 1 + 2a 



(5.12) 
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so that 

_ f 1 , , , 1 m 2 (l + 5a + 18a 2 + 20a 3 + 16a 4 ) / 

Im / s {x)dx = — + — + (m 

Jo 1 — a 2 8 1 — a 2 d 1 + 2a 1 + 3a v 



a 2 ) 3 (l + 2a)(l + 3a) 

(5.13) 

which yields the spectral density in the limit p —>■ oo. This result for the 
unperturbed term was first derived in [8]. 

We are now in a position to evaluate the perturbation term. The integral 
D can be computed from (I4.12p as 

1 m(2 + 4a + 5a 2 ) 



2(1 -a) 3 (l + a) 2 2(l-a) 4 (l + a) 3 (l + 2a) 
m 2 (5 + 35a + 128a 2 + 252a 3 + 256a 4 + 144a 5 ) ^ ( , 

+ 32 ; — — — r + O [m 6 

16(l-a) 5 (l + a) 4 (l + 2a) 2 (l + 3a) V 

(5.14) 

Then it follows from ( T4T71) that 

1 + 2a + 2a 2 



Im / si(x)dx 
Jo 



+ 



2(1 + 2a)(l - a 2 ) 3 
3m 2 (5 + 55a + 314a 2 + 1034a 3 + 2068a 4 + 2648a 5 + 1920a 6 + 768a 7 ) 
16(1 + 2a) 2 (l - a 2 ) 5 (l + 7a + 12a 2 ) 

+ 0(m 4 ). (5.15) 

As mentioned in the Introduction, Rodgers and BrayjTO] analysed sparse 
random matrices and derived a formula for the asymptotic spectral density 



pM~^{(m c ) 2 -^} 1/2 



(5.16) 



with /i 2 = 4{p + 1 + 0(l/p)}. Sparse random matrices in [10] give the 
adjacency matrices of Erdos and Renyi's random graphs, which correspond 
to the limit a — ► of Goh, Kahng and Kim's complex networks. From (15.131) 
and (I5.15p . we obtain 

rl 1 / \ 

Im / s (x)dx = -1 + -m 2 + O (m 4 ) (5.17) 
Jo 8 v ' 

Im [ sx(x)dx = — + — m 2 + (m 4 ) (5.18) 
Jo 2 16 v ' 



and 
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in the limit a — > 0. It follows that 

'<"> ~^ P [ 1 -T + M + \{\~ + ° m} + ° K 

(5.19) 

in agreement with (15.161) . as expected. 

6 Tail Region 

In order to analyse the EMA spectral density in the tail region, we need 
to work out the asymptotic formulae in the limit m — > oo for the integrals 
gUD and (Q5jl . 

To begin with, we define c = a — ib (with real a and 6) as 

c= (1 -a) 2 [ s (x)x- a dx (6.1) 
Jo 

and rewrite (14.4ft as 



s (x)x~ a = . (6.2) 

y ' mx a - c K ' 



Integrating over x from to 1, one obtains 

-I i 



-dx. (6.3) 



(1 — a) 2 Jo mx a — c 

This integral equation determines the parameter c. In [8], asymptotic for- 
mulae 

1 — a 



~ 7T- 



m 

(1 _ a )(l+«)/a J 



(6.4) 



a m (2-a)/a 

were derived in the limit m — > oo. It follows that 

Im / s (x)dx = Im / dx 

Jo Jo mx a — c 

(\-a) l/a 1 
-A— (6.5) 

Q m (2+a)/«' v I 

which gives the tail behaviour of the spectral density in the limit p — > oo. 
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Let us evaluate the 1/p correction to this asymptotic formula. Using the 
expressions 

J i = / jzzzr n — I^ dx > J 2 = / 7z— n — I^ dx > ( 6 - 6 ) 



o (mx a — c) 2 ' Jo (mx a — c) 



Jl = / 7 a ^3 dX ' J2 = / 7 a \3 dX 

Jo [mx a — c) 6 jo [mx a — c) a 

and (16.41) . we can derive the asymptotic formulae 

Re Ji ' ' 



1 — a m 2 ' 



(1 - a)^ 1 -^ 1 . . 

Im Ji ~ -7T 6.8 



Re Ji 



1 — a m 3 ' 



2 a 3 m (2+ a )/«' v ; 



Im J 2 ~ -tt^ ^- . . (6.10) 

a 2 m 2{1 - a )l a 



and 



7T (1 - a) (1_a)/Q 1 

Im J 2 ^ — — — 6.11 

2 a 3 m( 2 ~ Q )/ Q v ; 

in the limit m — > oo. Moreover, omitting logarithmic factors of the form 
mm, we find the asymptotic estimates 

Re/2 \ OK^V"), a>l/2 (6 - 12) 

and 

/ 0(m" 3 ), «<l/2, 

Re Js ~ { 0(m(«- 2 V«), a > 1/2. (6 - 13) 

Putting the above asymptotic formulae into (14.131) . we obtain 

f 1 , . , (l + «) (l-a) 1/a 1 

Im / si x dx ~ -tt^ ^ ^ p— . 6.14 

Jo « 3 m (2+3a)/a V 7 
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7 Discussion 

Using the replica method and effective medium approximation (EMA) in sta- 
tistical physics, we evaluated the spectral density for the adjacency matrix of 
Goh, Kahng and Kim's model of complex networks with a finite mean degree 
p. The EMA result was compared with numerically generated spectral den- 
sity. In the limit p —>■ oo, known results derived in [8] were reproduced. More- 
over, perturbative analytic formulae were presented for the EMA solution in 
the forms of 1/p expansions. As shown in Figure 1, although the agreement 
of the EMA spectral density with numerical data was fairly good, a signifi- 
cant discrepancy was observed around the band center. In order to improve 
the agreement, it seems necessary to develop more accurate schemes, such as 
a non-perturbative technique |10j or single defect approximation [Til [T2| [14"]. 
The discrepancy around the band center might be related to the occurrence 
of the eigenvector localisation. In connection to the localisation, we expect 
that further studies on the local properties of the spectra, such as the dis- 
tribution of the eigenvalue spacings[H)| or the largest eigenvalues [TB"]. will be 
illuminating. 
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